#!/usr/bin/perl
use Cwd 'abs_path';
#use Getopt::Std;
#Input Arguments
#getopt("d:t:a:f:o", \%args);
#my $dir_or = $args{d};
use Getopt::Long;
&Getopt::Long::GetOptions(
'config=s'      => \$config
);
$line = abs_path($0);
chomp $line;

@DR_array = split('/',$line);
pop(@DR_array);
$dir = join("/",@DR_array);
require "$dir/bin/CONFIG.pl";
getDetails($config);
#my $tped= $args{t};
#my $arraytype= $args{a};
#my $tfam = $args{f};
#my $outdir = $args{o};
#my $dirtemp = $args{d};
my $tped= $config{"TPED"};
my $tfam = $config{"TFAM"};
my $outdir = $config{"OUTPUT_FOLDER"};
my $dirtemp = $config{"TEMP_FOLDER"};
my $forward_ind = $config{"FORWARDSTRAND_IND"};
print "TPED:$tped\nTFAM:$tfam\nOUTDIR:$outdir\nDIRTEMP:$dirtemp\nFORWAD STRAND IND: $forward_ind\n";
chomp($dirtemp);
chomp($tped);
chomp($tfam);
chomp($outdir);
if($dirtemp eq "")
{
        die "entered temp dir name  is empty\n";

}
if (!(-e $tfam)) {
        die "input tfam file does not exist\n";
}
if (!(-e $tped)) {
        die "input tped file does not exist\n";
}
if (!(-e $forward_ind)) {
        die "Forward strand IND does not exist\n";
}

if($tped eq "")
{
        die "entered tped file name  is empty\n";

}
else
{
@tped_ar=split(/\//,$tped);
$tped = pop(@tped_ar);
$dir_or = join('/',@tped_ar);
}
if($tfam eq "")
{
	die "entered tfam file name  is empty\n";
}
if($outdir eq "")
{
	die "entered output directory name  is empty\n";
} 

#to know the path of this script to grab working directory location
#$line = abs_path($0);
#chomp $line;

#@DR_array = split('/',$line);
#pop(@DR_array);
#$dir = join("/",@DR_array);
$round = sprintf("%.0f", rand()*time());
$rounded = "temp".$round;
#$rounded = "temp545278769";
#creating a temp directory if not presen
system("mkdir -p $dirtemp");
system("mkdir $dirtemp/$rounded");
system("cp $tfam  $dirtemp/$rounded");
@tfam_ar=split(/\//,$tfam);
$tfam = pop(@tfam_ar);
system("mv $dirtemp/$rounded/$tfam $dirtemp/$rounded/unprocessed_beagle_input.tfam");
system("cp $dir_or/$tped $dirtemp/$rounded/unprocessed_beagle_input.tped");
$tfam ="unprocessed_beagle_input.tfam";
$tped = "unprocessed_beagle_input.tped";

system("cp $dirtemp/$rounded/unprocessed_beagle_input.tfam $dirtemp/$rounded/processed_beagle_input.tfam");
$command_sys = "/usr/bin/perl $dir/bin/process_fwdstnd_beagle_input.pl -f $dirtemp/$rounded/unprocessed_beagle_input.tped -h $dir/bin/files/newhapmap3_sample_run.tped  -n $dirtemp/$rounded/processed_beagle_input.tped -e $dirtemp/$rounded/excluded_no_hapmap_processed_beagle_input.tped -a $dirtemp/$rounded/ambi_hapmap_processed_beagle_input.tped -u $dirtemp/$rounded/unsure_differentallele.tped  -i $forward_ind -r $dir/bin/files/snps_notin_xchr_beagle_imp_reference -o $dirtemp/$rounded/snps_notin_xchr_beagle_imp_reference_output";
print $command_sys."\n";
system($command_sys);
$command_sys = "/usr/bin/perl $dir/bin/perl_snps_notref_dosage.pl $dirtemp/$rounded/excluded_no_hapmap_processed_beagle_input.tped  $dirtemp/$rounded/snps_noref_tped.dosage";
print $command_sys."\n";
system($command_sys);
system("$dir/bin/plink/plink --tfile $dirtemp/$rounded/processed_beagle_input --out $dirtemp/$rounded/processed_beagle_input --make-bed");
#chopping the chr 1-22 and submitting the jobs for imputation
$syst = "/usr/bin/perl $dir/bin/perl_create_ind_chr_imputation_ori.pl -d $dir -t $dirtemp/$rounded";

print $syst."\n";
system($syst); 
$sys = "/usr/bin/perl $dir/bin/perl_create_ind_chr_imputation_ori_X.pl -d $dir -t $dirtemp/$rounded";
#submitting the chr x
system($sys);
open(BU,"$dirtemp/$rounded/submitjobs");
open(WRBU,"$dirtemp/$rounded/submitjobs1");
$line_jobs = "";
$submitjobs_num = 0;
while($line = <BU>)
{
	chomp($line);
	$submitjobs_num++;
	$line = $line." >>$dirtemp/$rounded/submitjobs1";
	print $line."\n";
	system($line);
	
}
$submitjobs_num1 = 0;
open(BU1,"$dirtemp/$rounded/submitjobs1");
while($line = <BU1>)
{
        chomp($line);
	$submitjobs_num1++;
	@array = split(" ",$line);
	$line_jobs = $line_jobs.$array[2]." ";
}
if($submitjobs_num != $submitjobs_num1)
{
	die "all jobs are not submitted\n";
}

#checking the submitted jobs which are completed
system("JAVA_HOME=/usr/java/latest");
system("export JAVA_HOME");
system("/home/inforsense/applications/groovy-1.6.4/bin/groovy  $dir/bin/CheckJobStatus.groovy $line_jobs");
#merging the result files chromosome wise 
#open(WWW, ">$dic/temp/$rounded/failed.txt");
open(WWW, ">$dirtemp/$rounded/failed.txt");
$sys= "/usr/bin/perl $dir/bin/perl_after_imputation_merge.pl -d $dir -t $dirtemp/$rounded";
print $sys."\n";
system($sys);
system("wait");
$sys = "/usr/bin/perl $dir/bin/perl_xchr_convert_tped_create_dosage.pl -d $dirtemp/$rounded";
print $sys."\n";
system($sys);

#creating the combined male and female  x chr imputation results files
###system("/usr/bin/perl $dir/bin/perl_club_male_female_x_chr_after_imp.pl -d $dir/temp/$rounded");
open BUF,"< $dirtemp/$rounded/failed.txt" or die "can't open failed.txt";
open(BBUFFB,">$dirtemp/$rounded/beagle_errors");
close(BBUFFB);
#checking if any chr outfile is not created or failed
while($lined = <BUF>)
{
	if($lined =~ m/failed/)
	{
		system("rm $dirtemp/$rounded/whole_geno_phased");
		system("rm $dirtemp/$rounded/whole_gprobs");
		system("rm $dirtemp/$rounded/whole_no_snp_hapmap.txt");
		system("rm $dirtemp/$rounded/whole_r2");
		@kgk = split(/\//,$lined);
		pop(@kgk);
		$linedd = join('/',@kgk);
		system("cat $linedd/run*.e* >> $dirtemp/$rounded/beagle_errors");
		die "$lined\n";
	}
}
$maindir = "$dirtemp/$rounded";
system("/usr/bin/perl $dir/bin/perl_create_beagle_output_plink_tped.pl -m $dir/bin/files/newhapmap3_sample_run.tped -i $maindir/whole_geno_phased -o $maindir/whole_geno_phased.tped");
#system("/usr/bin/perl $dir/bin/perl_create_combine_dosage_1stphase_nox.pl $maindir/whole_gprobs $maindir/whole_gprobs1 $maindir/whole_geno_phased.tped");
$maindir = "$dirtemp/$rounded";
open(ALLF,">$maindir/allfiles.txt");
$temp_kk =0;
if(-e "$maindir/23/male/geno_male.bgl.gz") 
{
	$temp_kk =1;
}
if(-e "$maindir/23/male/geno_male.bgl")
{
	$temp_kk =1;
}
print "male $temp_kk $maindir/23/male/geno_male.bgl.gz\n";
if($temp_kk ==1)
{

        $file = "$dirtemp/$rounded/23/male/male.tped";
        if (!(-e $file)) 
	{
                die "male tped file does not exist\n";
        }
        system("$dir/bin/plink/plink --tfile $dirtemp/$rounded/23/male/male --make-bed --out $dirtemp/$rounded/23/male/male");
	print ALLF "$maindir/23/male/male.bed $maindir/23/male/male.bim $maindir/23/male/male.fam\n";
	#system("gunzip $maindir/23/male/*");
	system("wc -l $maindir/23/male/male.tped > $maindir/dosage_x.size");
}

if(-e "$maindir/23/female/geno_female.bgl.gz")
{
        $temp_kk =2;
}
if(-e "$maindir/23/female/geno_female.bgl")
{
	$temp_kk =2;
}
print "female $temp_kk $maindir/23/female/geno_female.bgl.gz\n";
if($temp_kk ==2)
{


        $file = "$dirtemp/$rounded/23/female/female.tped";
        if (!(-e $file)) 
	{
                die "female tped file does not exist\n";
        }
        system("$dir/bin/plink/plink --tfile $dirtemp/$rounded/23/female/female --make-bed --out $dirtemp/$rounded/23/female/female");
	print ALLF "$maindir/23/female/female.bed $maindir/23/female/female.bim $maindir/23/female/female.fam\n";
	#system("gunzip $maindir/23/female/*");
	system("wc -l $maindir/23/female/female.tped > $maindir/dosage_x.size");
}
close(ALLF);
system("cp $maindir/processed_beagle_input.fam $maindir/whole_geno_phased.tfam");
if($temp_kk !=0)
{
	system("$dir/bin/plink/plink --tfile $maindir/whole_geno_phased --merge-list $maindir/allfiles.txt --transpose --recode --out $maindir/whole_geno_phased_x");
	open(SIZE,"$maindir/dosage_x.size");
	$dosage_size = <SIZE>;
	@dosage_size = split(" ",$dosage_size);
	$dosage_size =  $dosage_size[0];
	$sys = "tail -$dosage_size $maindir/whole_geno_phased_x.tped > $maindir/dosage.tped";
	print $sys."\n";
	system($sys);
	if($temp_kk == 2)
	{
		$sys = "/usr/bin/perl $dir/bin/perl_xchr_create_dosage.pl -d $dirtemp/$rounded";
	}
	else
	{
		$sys = "/usr/bin/perl $dir/bin/perl_xchr_create_dosage_nofemale.pl -d $dirtemp/$rounded";
	}
	print $sys."\n";
	system($sys);
	$maindir = "$dirtemp/$rounded";
	system("cat $maindir/whole_gprobs $maindir/X_dosage> $maindir/whole_gprobs_x");
	system("cp $maindir/whole_geno_phased_x.tped $maindir/whole_geno_phased_x_ori.tped");
	system("cp $maindir/whole_gprobs_x $maindir/whole_gprobs_x_ori");
	system("mv $maindir/whole_geno_phased_x.tped $maindir/whole_geno_phased.tped");
	system("mv $maindir/whole_gprobs_x $maindir/whole_gprobs");
	$maindir = "$dirtemp/$rounded";
}
$maindir = "$dirtemp/$rounded";
system("/usr/bin/perl $dir/bin/perl_create_combine_dosage_1stphase_nox.pl $maindir/whole_gprobs $maindir/whole_gprobs1 $maindir/whole_geno_phased.tped");
system("wc -l $maindir/ambi_hapmap_processed_beagle_input.tped > $maindir/ambi_count");
open AMBI,"$maindir/ambi_count" or die "nor file found $maindir/ambi_count";
$line = <AMBI>;
@arr = split(" ",$line);
print "number of lines in ambi file $arr[0]\n";
if($arr[0] > 0)
{
system("/usr/bin/perl $dir/bin/perl_check_ambi_snp.pl -a $maindir/ambi_hapmap_processed_beagle_input.tped -i $maindir/whole_geno_phased.tped -o $maindir/whole_geno_phased_out.tped -d $maindir/whole_gprobs1 -n $maindir/whole_chr_gprobs_dosage_ambi");
system("mv $maindir/whole_geno_phased_out.tped $maindir/whole_geno_phased.tped");
system("mv $maindir/whole_chr_gprobs_dosage_ambi $maindir/whole_gprobs1");
}
close(AMBI);
system("/usr/bin/perl $dir/bin/perl_create_indicator_file_1stphase_nox.pl $maindir $maindir/whole_r2 $maindir/unprocessed_beagle_input.tped $maindir/whole_geno_phased.tped $maindir/whole_gprobs1");
$maindir = "$dirtemp/$rounded";
system("/usr/bin/perl $dir/bin/perl_add_snps_not_in_impref_2results.pl $maindir/snps_notin_xchr_beagle_imp_reference_output $maindir/whole_gprobs1 $maindir/whole_ind $maindir/whole_geno_phased.tped");
#$syste = "/usr/bin/perl $dir/bin/perl_create_new_map_tpedfile_4sort.pl $maindir/whole_geno_phased.tped $maindir/snp_order.txt $maindir";
#print $syste."\n";
#system($syste);
#system("/usr/bin/perl $dir/bin/perl_inputlist_sortfile.pl $maindir/snp_order.txt $maindir/whole_gprobs1 1 $maindir/sorted_whole_gprobs1");
#system("/usr/bin/perl $dir/bin/perl_inputlist_sortfile.pl $maindir/snp_order.txt $maindir/whole_ind  1 $maindir/sorted_whole_ind");
#system("/usr/bin/perl $dir/bin/perl_inputlist_sortfile_tped.pl $maindir/snp_order.txt  $maindir/whole_geno_phased.tped 2 $maindir/sorted_whole_geno.map");
$maindir = "$dirtemp/$rounded";
system("cut -f1-4 -d ' ' $maindir/whole_geno_phased.tped|perl -e 'while(<STDIN>){\$_ =~ s/^23 /X /g;\$_ =~ s/^24 /Y /g;\$_ =~ s/^26 /M /g;print \$_}'> $maindir/extracted_whole_geno.map");
system("sort -k1,1 -k4,4n -k2,2 -T $maindir $maindir/extracted_whole_geno.map > $maindir/sorted_whole_geno.map");
system("rm $maindir/extracted_whole_geno.map");
$syste = "/usr/bin/perl $dir/bin/perl_final_4sort.pl $maindir/sorted_whole_geno.map $maindir/whole_gprobs1 $maindir/sorted_whole_gprobs1 $maindir/whole_ind  $maindir/sorted_whole_ind";
print $syste."\n";
system($syste);

#system("$dir/bin/sort -k1,1 -k4,4n -k2,2 $maindir/whole_geno_phased.tped |cut -f1-4 -d ' ' |perl -e 'while(<STDIN>){\$_ =~ s/^23 /X /g;\$_ =~ s/^24 /Y /g;\$_ =~ s/^26 /M /g;print \$_}'> $maindir/sorted_whole_geno.map");
system("sort -k1,1 -k4,4n -k2,2 -T $maindir $maindir/sorted_whole_gprobs1 |sed 's/^[A-Za-z0-9]* [a-z0-9]* [a-z0-9]* [a-z0-9]* //'> $maindir/sorted_whole_gprobs2");
system("mv $maindir/sorted_whole_gprobs2 $maindir/sorted_whole_gprobs1");

system("sort -k1,1 -k4,4n -k2,2 -T $maindir $maindir/sorted_whole_ind |sed 's/^[A-Za-z0-9]* [a-z0-9]* [a-z0-9]* [a-z0-9]* //'> $maindir/sorted_whole_ind1");
system("mv $maindir/sorted_whole_ind1 $maindir/sorted_whole_ind");
if (-e "$outdir") 
{
	print "Output directory already exists: $outdir";
}
else 
{
	system("mkdir -p $outdir");
}
system("sort -k1,1 -k4,4n -k2,2 $maindir/snps_noref_tped.dosage > $maindir/snps_noref_tped1.dosage");
#code to remove snps in the same location

open(BUFF,"$maindir/snps_noref_tped1.dosage") or die " no file found $maindir/snps_noref_tped1.dosage\n";
open(WRBUFF,">$maindir/snps_noref_tped.dosage") or die " not able to write $maindir/snps_noref_tped.dosage\n";
open(WRBUFF1,">$maindir/final_removed.extra") or die "not able to write $maindir/final_removed.extra\n";
@array = ();
$prevsnp = "";
while(<BUFF>)
{
        chomp($_);
        @a = split(" ",$_);
        if($prevsnp ne "" && $prevsnp ne "$a[0] $a[3]" )
        {
                #print @array."\n";
                $num = 0;
                for($i=0;$i<@array;$i++)
                {
                        @b = split(" ",$array[$i]);
                        if($num == 0 &&($b[1] =~ m/^rs/ || $i == @array-1))
                        {
                                print WRBUFF $array[$i]."\n";
                                $num++;
                        }
                        else
                        {
                                print WRBUFF1 $array[$i]."\n";
                        }
                }
                @array = ();
        }
        unshift(@array,$_);
        $prevsnp = "$a[0] $a[3]";
}
$num = 0;
for($i=0;$i<@array;$i++)
{
        @b = split(" ",$array[$i]);
        if($num == 0 &&($b[1] =~ m/^rs/ || $i == @array-1))
        {
                print WRBUFF $array[$i]."\n";
                $num++;
        }
        else
        {
                print WRBUFF1 $array[$i]."\n";
        }
}
system("sort -k1,1 -k4,4n -k2,2 $maindir/snps_noref_tped.dosage > $maindir/snps_noref_tped1.dosage");
system("mv $maindir/snps_noref_tped1.dosage $outdir/final.extra");
#system("mv $maindir/snps_noref_tped1.dosage $maindir/snps_noref_tped.dosage");
system("mv $maindir/sorted_whole_gprobs1 $outdir/final.dosage");
system("mv $maindir/sorted_whole_ind $outdir/final.ind");
system("mv $maindir/sorted_whole_geno.map  $outdir/final.map");
if(-e "$maindir/whole_geno_phased_x.tfam $outdir/final.tfam")
{
	system("mv $maindir/whole_geno_phased_x.tfam $outdir/final.tfam");
}
else
{
	system("mv $maindir/whole_geno_phased.tfam $outdir/final.tfam");
}
	#system("mv $maindir/final_removed.extra $outdir/final_removed.extra");
#system("mv $maindir/samples.qc $outdir/sample.qc");
#system("rm -rf $maindir/*");
#system("rm -rf $maindir");
